	subroutine meshCellOut_def_old(xg,n,meshCellOut,nCellsOut,i,tol)
	implicit none
	real *8 xg(3)
	integer nCellsOut,i
	real *8 n,tol
	integer meshCellOut(*)
	real *8 d,rr1,rr2,vtmp
	vtmp=0.14d0*n
	rr1=vtmp-tol
	rr2=vtmp+tol
	if (xg(1).gt.rr1.and.xg(1).lt.rr2 ) then
	if (dabs(xg(2)).lt.10d-4.or.dabs(xg(3)).lt.10d-4) then
        nCellsOut=nCellsOut+1
        meshCellOut(nCellsOut)=i
	return
        end if
	if (dabs(xg(3)).lt.10d-4.or.dabs(xg(2)).lt.10d-4) then
        nCellsOut=nCellsOut+1
        meshCellOut(nCellsOut)=i
        end if
        end if
	return
        end 

	subroutine meshCellOut_def(xg,n,meshCellOut,nCellsOut,i,tol)
	implicit none
	real *8 xg(3)
	integer nCellsOut,i
	real *8 n,tol,dtol
	integer meshCellOut(*)
	real *8 d,rr1,rr2,vtmp
	vtmp=0.14d0*n
	dtol=5d-3
	rr1=vtmp-tol
	rr2=vtmp+tol
	if (xg(1).gt.rr1.and.xg(1).lt.rr2 ) then
	if (dabs(xg(2)).lt.dtol.or.dabs(xg(3)).lt.dtol) then
        nCellsOut=nCellsOut+1
        meshCellOut(nCellsOut)=i
	return
        end if
	if (dabs(xg(3)).lt.dtol.or.dabs(xg(2)).lt.dtol) then
        nCellsOut=nCellsOut+1
        meshCellOut(nCellsOut)=i
        end if
        end if
	return
        end 
	subroutine meshCellOut_defn(x,meshCellOut,nCellsOut,iel, &
	area,xCenter,meshSf)
	implicit none
	real *8 area,xCenter(3),meshSf(3)
	integer nCellsOut,iel
	integer meshCellOut(*)
	integer rib(2,12)
	real *8 x(3,8)
	real *8 xrib(3,12),xCells(3,12)
	integer location,i,j,k,locationp,nPoints,nVar(4,6),n1,n2,n3,n4
	real *8 intersect(12),mag,z(3)
	data (rib(1,i),i=1,12)/1,2,3,4,1,2,3,4,5,6,7,8/
	data (rib(2,i),i=1,12)/2,3,4,1,5,6,7,8,6,7,8,5/
	data (nVar(1,i),i=1,6)/1,1,1,1,1,1/
	data (nVar(2,i),i=1,6)/2,2,3,3,4,4/
	data (nVar(3,i),i=1,6)/3,4,2,4,2,3/
	data (nVar(4,i),i=1,6)/4,3,4,2,3,2/
!	write(*,*) x
	call surface(x(1:3,1),location)
	do i=2,8
	locationp=location
	call surface(x(1:3,i),location)
	if (locationp.ne.location) then
!	write(*,*) nCellsOut+1,iel
        nCellsOut=nCellsOut+1
        meshCellOut(nCellsOut)=iel
	nPoints=0
	xCenter=0d0
	do j=1,12
	call line(x(1:3,rib(1,j)),x(1:3,rib(2,j)),xrib(1:3,j),intersect(j))
	if (intersect(j).ge.0d0) then
	nPoints=nPoints+1
	xCells(1:3,nPoints)=xrib(1:3,j)
	xCenter(1:3)=xCenter(1:3)+xrib(1:3,j)
	write(*,*) j,intersect(j)
	end if
	end do
	xCenter=xCenter/nPoints

    call sxdef(meshSf(1:3),xCells(1:3,1),xCells(1:3,2),xCenter(1:3),xCenter(1:3))
	mag=dsqrt(meshSf(1)**2+meshSf(2)**2+meshSf(3)**2)
	meshSf=meshSf/mag
	z=xCenter+meshSf*1d-3
	call surface(z,location)
	if (location.lt.0) meshSf=-meshSf
	if (nPoints.gt.4) stop 'Fatal Error in FW'
	if (nPoints.eq.3) then
    call sxdef(z,xCells(1:3,1),xCells(1:3,2),xCells(1:3,3),xCells(1:3,3))
	area=dsqrt(z(1)**2+z(2)**2+z(3)**2)
	return
	end if
	area=0d0
	do k=1,6
	n1=nVar(1,k)
	n2=nVar(2,k)
	n3=nVar(3,k)
	n4=nVar(4,k)
    call sxdef(z,xCells(1:3,n1),xCells(1:3,n2),xCells(1:3,n3),xCells(1:3,n4))
	mag=dsqrt(z(1)**2+z(2)**2+z(3)**2)
	if (area.lt.mag) area=mag
	end do
	


!	stop
	return
	end if
	end do

	return
	end
	subroutine surface(x,n)
	implicit none
	real *8 x(3)
	real *8 rr,rr1
	integer n
	n=1
	rr=dsqrt(x(2)**2+x(3)**2)
	rr1=0.17 !*28
!	rr1=0.158+0.158*14.0/180.0*3.1415927*x(1)
	if (rr.lt.rr1) n=-1
	return
	end
	subroutine line(x,y,z,alpha)
	implicit none
	real *8 x(3),y(3),z(3),alpha
	real *8 x0,x1,x2
	integer loc,loc1,loc2
	x1=0d0
	x2=1d0
	call surface(x,loc1)
	call surface(y,loc2)
	alpha=-1d0
	if (loc1.eq.loc2) return
 1	x0=(x1+x2)/2d0
 	z(1:3)=x(1:3)*(1.0d0-x0)+y(1:3)*x0
	call surface(z,loc)
	if (loc.eq.loc1) x1=x0
	if (loc.eq.loc2) x2=x0
	if (dabs(x1-x2).gt.1d-6) goto 1
	alpha=x0
	
	return
	end

